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Chapter  1 


Imaging  in  a  parallel-plate  waveguide 


This  is  joint  work  with  Cliff  Nolan.  Work  on  this  project  is  still  ongoing. 


1.1  Basic  theory 


The  mathematical  model  involves  a  number  of  ingredients:  1)  a  model  for  wave  propagation  in 
free  space,  2)  a  model  for  multiple  scattering  from  buildings,  3)  a  (linearized)  model  for  scattering 
from  the  scene,  and  4)  a  model  for  the  array  of  sources. 


Figure  1.1:  Scattering  geometry.  The  darker  walls  and  dot  denote  the  true  walls  and  a  point  source; 
the  lighter  ones  represent  virtual  walls  and  corresponding  virtual  sources. 


1.1.1  Model  for  wave  propagation  in  free  space 

We  assume  that  waves  propagate  within  the  waveguide  according  to  the  scalar  wave  equation: 

V2u  —  Cq  2ii  =  0,  (1-1) 

where  the  dots  denote  differentiation  with  respect  to  t  and  c0  is  the  speed  of  light. 
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(1.2) 


The  eld  due  to  a  point  source  at  P  =  0  at  time  t  —  0  is 


9o(t,\P\)  = 


S(t-\P\/c0 ) 
4tt|P| 


3-iw(t-|P|/co) 


4tt|P| 


-dcu 


We  will  use  capital  letters  for  frequency-domain  quantities,  which  are  related  to  time-domain 
quantities  by  the  Fourier  transform 

U (u,  P)  =  —  [  eiu>tu(t,  P)dt.  (1.3) 

27T  J 

We  write  k  —  uj/cq. 


1.1.2  Wave  propagation  in  the  waveguide 


We  consider  the  waveguide  formed  by  two  in  nite  parallel  walls.  We  take  one  wall  to  be  the  plane 
x  =  0  and  the  other  to  be  the  plane  x  =  L.  We  assume  that  on  the  walls,  the  electromagnetic  eld 
satisi  es  the  boundary  condition 


^  =  0 
OX 


(1.4) 


Consequently  a  Green’s  function  for  the  waveguide  can  be  constructed  by  the  method  of  images: 


G(u>,P,P) 


£  £ 

ae{±) n=—o o 


e'k\P-K,a\ 
4vr|  P-P;, 


(1.5) 


where  the  P'n  ±  are  the  virtual  images  of  P'\  if  P'  =  ( x y\  z'),  then  P'l  +  =  ( x '  +  2 nL,  y' ,  z ')  and 
P'n  _  =  (— x'  +  2 nL,  y',  z')  (See  Fig.  3.1).  The  sum  in  (1.5)  is  over  all  the  virtual  images. 


We  note  that  the  virtual  images  satisfy  certain  symmetry  relations: 


1  p  —  p1  1  — 
l-1^  rn,+ 1 

y/[x  ~ 

-  ( x ’  +  2 nL)}2  +  •  •  • 

=  \P-n,+  -  P'l 

IP-P'  1  = 

1  n,—  1 

Vlx  - 

(—x1  +  2  nL)]2  +  •  • 

= 

Vi~x 

+  2nL  —  x')}2  +  •  ■  ■ 

=  lpn,_-p,l 

and  that 

(P  +  p)n,a  —  Pn,a  +  P0,a-  (1-7) 

We  use  tildes  to  indicate  re  ection  of  the  ^-coordinate  about  the  y-z- plane:  p  =  p0  _  and  = 

P  =  Po.+- 


Within  the  region  —  L  <  x  <  L,  G  satis  es 

V2G(co,  P,  P')  +  k2G(u,  P,  P’)  =  —5(P,  P’) 


(1.8) 


To  model  a  nite  waveguide,  say  one  extending  up  to  height  z  =  H,  we  include  only  certain  of 
the  virtual  sources. 
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1.1.3  The  model  for  scattering  from  the  target 


We  model  the  target  as  a  change  in  wave  speed  V(P)  =  c~2(P )  —  Cq2.  The  scattered  eld 
Usc(u,  P,  P')  at  P  due  to  an  incident  eld  Um(u,  P,  P')  =  G(u,P,P ')  is  given  by  [31]  the 
Lippmann-Schwinger  equation: 

Usc(tu,P,P')  =  I  G(u,P,P")cu2V(P")[G(u,P",P')  +  Usc(u,P",P')}dP"  (1.9) 
We  wish  to  reconstruct  V  from  knowledge  of  Gsc. 

Unfortunately  ( 1 .9)  is  an  integral  equation  whose  solution  Usc  depends  in  a  nonlinear  way  on 
V.  To  simplify  the  problem,  we  make  the  Born  or  single -scattering  approximation:  we  neglect  the 
term  involving  Usc  on  the  right  side  of  (1.9).  This  gives  us 

Ug(u,P,P')=  I  G(u,  P,  P" )u2V{P")G(u,  P",  P’)dP"  (1.10) 


The  Born  approximation  makes  the  mapping  from  V  to  Usc  linear,  but  it  is  not  necessarily 
a  good  approximation.  Another  linearizing  approximation  that  can  be  used  for  re  ection  from 
smooth  surfaces  is  the  Kirchhojf  approximation,  in  which  the  scattered  eld  is  replaced  by  its 
geometrical  optics  approximation.  Here,  however,  we  consider  only  the  Born  approximation. 

The  corresponding  time-domain  eld  is 

us£(t,P,P')  =  [  e~'lultG{uj:P:P")V{P")G{uj,P",P')u2dudP"  (1.11) 


We  note  that  since  G  is  given  by  a  sum  ( 1 .5),  ( 1 . 1 1)  is  of  the  form 


us£(t,P,P')  =  Y.F^mAV](^P>P') 

a,/3e{±}  m,n 


(1.12) 


where 


r  e-iu,(HP-p"a|/co)  J“\p"~pLJ/c° 

Fn^mAV}(t,P,P')=  I  „„  V{p").  ^dwdP" 


4vr|P  -  P"c 


4tt|  P"-Pn 


m,(3 1 


(1.13) 


To  the  initial  factors  in  (1.13)  we  apply  the  symmetry  relations  (1.6)  and  then  we  re-label  the 
indices  of  ( 1 . 12)  to  write  usc  in  the  same  fonn  (1.12)  with 


F1wnAV}(t,P,P')  = 


^-iw{t-\Pn,a-P"\/co) 

47r|P„ia  -  P 


eMP"-PL,p\/c0 

-V(P")  — — - o2dudP" 


4tt|  P" 


r»  | 

rm,j3  I 


(1.14) 
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1.1.4  Arrays  of  sources 


We  consider  arrays  of  isotropically-radiating  sources.  At  frequency  u>,  the  eld  at  P"  emanating 
from  the  source  at  P'  is  G( u,  \P"  —  P'\),  where  G  is  given  by  (1.5).  To  model  the  eld  um  from 
an  array,  centered  at  P1,  in  which  the  element  at  P'  +  p'  is  activated  with  the  waveform  j(t,p'), 
we  simply  integrate: 

/OO  p 

e~iuJt  /  G(cu,  | P"  -  ( P 1  +  pf)  | )  J {oj i  p'^dp1  du.  (1.15) 

-oo  J  array 

(Here  J  denotes  the  Fourier  transform  of  j.) 

To  obtain  the  corresponding  scattered  eld,  we  apply  the  same  operation  to  usc  (because  the 
mapping  from  incident  eld  to  scattered  eld  is  linear).  We  denote  the  resulting  Born-approxiimated 
eld  by  u“B: 

«“b  (t,P,P)  =  J  e~MG(u>,P,P")V(P") 

f  G(u>,  P",  P'  +  p')  J{oj,p')dp'oj2d(jj  (1.16) 

J  array 


An  array  can  form  a  steered  beam  by  including  phasing  in  j.  For  example,  for  a  unit  vector  p, 
abeam  steered  in  direction  p  and  be  produced  by  taking  p')  =  J(u>)  exp  ( i  k  p,  ■  p' )  ;vari.ay  ip' ) , 
where  xarray  denotes  a  function  that  is  smooth  on  the  array  and  zero  off  it. 


We  assume  that  p'  is  small  compared  to  | P"  —  P'\,  which  enables  us  to  use  the  expansion 

\P"  -  {P1  +p)\  =  \P"  -P'\  p  +  •••  (1.17) 


in  the  Fourier- transformed  version  of  (1.15): 


u'-(u,p",p')=  V  v 

ae{±}  n 
"  e 


jk(\p"-p^a-P'0j 


EE 

a£{±}  n  ' 


4tt|  P"  -  -  p',ia: 


An  |  P"  —  P'n 


—  J{u)e-^x^y{p')dp' 

J(ou)e-ik»p,xaTray(p')dp'  (1.18) 


Arrays  are  often  planar;  in  this  case  the  p'  integration  is  two-dimensional.  For  a  planar  array  and 
a  =  +,  the  integral 


3-i  KP"-Pk,a-Pro,+-^P/)dpl  =  /  e-i  dp’ 


(1.19) 


'  array 


'  array 


appearing  in  (1.18)  is  a  product  of  sine  functions  whose  main  lobe  is  at  P"  —  P'na  =  p.  When 
a  =  — ,  the  same  integral 


l-\k{P"-p^a-P/0_-^v')dpt  =  /  e_i  k{P"-p^a--H)-pf  dpt 


(1.20) 


'  array 


1  array 
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has  its  main  lobe  in  direction  P"  —  P'na  =  [i 

In  some  of  the  cases  to  follow,  it  is  important  to  be  able  to  form  a  beam;  in  others  it  is  not.  In 
the  cases  in  which  it  is  not,  we  consider  only  the  eld  emanating  from  a  single  source. 

Reception  by  the  array  is  modeled  by  integrating  the  scattered  eld  over  the  array,  again  perhaps 
with  some  weight  for  steering. 


1.2  Monostatic  imaging 


Work  on  this  part  of  the  project  is  ongoing,  but  is  close  to  being  nish  ed.  The  conclusion  is  that 
the  standard  monostatic  SAR  scenario  will  give  rise  to  images  that  have  an  unacceptable  number 
of  artifacts.  However,  a  more  workable  scenario  is  to  use  a  x  ed  array  antenna  that  produces  a 
beam  suf  ciently  narrow  to  distinguish  among  different  scattering  paths.  Such  measurements  can 
be  used  to  produce  images  free  from  artifacts. 


1.3  Bistatic  imaging 


We  have  considered  two  scenarios:  a)  a  x  ed  transmitter  (perhaps  a  source  of  opportunity)  and  a 
moving  receiver,  and  b)  two  independently  moving  antennas  that  can  both  transmit  and  receive. 
Our  conclusions  are  that  both  scenarios  lead  to  images  with  an  unacceptable  number  of  artifacts. 
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Chapter  2 


Imaging  from  high-Doppler-resolution 
waveforms 


This  is  joint  work  with  Brett  Borden.  It  has  been  published  in  [6]. 


2.1  Introduction 


Standard  synthetic-aperture  radar  imaging  systems  [42]  transmit  wideband  waveforms,  and  the 
corresponding  radar  return  signals  are  processed  to  synthesize  the  response  from  a  sharp  delta-like 
pulse.  Such  wideband  waveforms  are  called  high-range-resolution  waveforms  because  their  radar 
returns  can  be  used  to  obtain  accurate  estimates  of  the  distance  (range)  to  a  scatterer. 


Figure  2. 1 :  A  high-range-resolution  system. 
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When  a  high-range-resolution  system  is  used  to  image  the  scatterers  on  a  at  surface,  the  radar 
return  at  each  time  t  is  a  superposition  of  all  the  returns  due  to  those  scatterers  positioned  at 
distance  2 t/c  from  the  radar.  The  imaging  problem  can  then  be  formulated  [22,  1,3,  26]  in  terms 
of  reconstructing  the  scattering  density  function  p  from  its  integrals  over  all  circles  centered  on  the 
ight  path  of  the  antenna  (see  Figure  3.1). 

Radar  systems  can  be  designed,  however,  to  operate  in  a  complementary  mode:  instead  of 
transmitting  high-range-resolution  pulses  and  estimating  target  range,  they  can  transmit  a  high- 
Doppler-resolution  waveform  (a  x  ed-frequency  waveform,  also  called  a  continuous-wave  or  CW 
waveform)  and  estimate  the  relative  target  velocity  from  the  Doppler  frequency  shift  of  the  return. 
For  an  antenna  moving  at  a  constant  velocity  over  a  at  surface,  the  return  at  a  given  Doppler  shift 
is  a  superposition  of  all  the  returns  due  to  scatterers  with  the  same  relative  velocity,  and  all  the 
scatterers  with  the  same  relative  velocity  lie  on  a  certain  hyperbola,  called  an  iso-Doppler  curve 
or  isodop.  This  suggests  a  corresponding  imaging  approach:  reconstruct  the  scattering  density 
function  p  from  its  integrals  over  the  iso-Doppler  hyperbolas  (see  Figure  3.2). 


Figure  2.2:  A  high-Doppler-resolution  system. 

High-Doppler-resolution  imaging  systems  (which  we  refer  to  as  simply  Doppler  imaging  sys¬ 
tems)  require  only  a  relatively  simple  (inexpensive)  transmitter  and  may  thus  have  advantages  over 
high-range-resolution  systems  in  some  situations.  In  addition,  Doppler  imaging  systems  may  be 
useful  in  scenarios  in  which  the  radar  signals  must  penetrate  through  a  medium  with  a  frequency- 
dependent  attenuation. 

The  concept  of  Doppler  imaging  is  not  entirely  new.  The  notion  appears  in  [28] — in  the  con¬ 
text  of  a  rotating  target  and  stationary  radar — to  motivate  the  development  of  range-Doppler  or 
inverse-synthetic-aperture  imaging.  The  idea  is  also  implicit  in  the  theory  of  space-time  adaptive 
processing  [24].  Moreover,  reconstruction  of  a  scattering  density  function  from  its  integrals  over 
hyperbolas  was  developed,  for  wideband  range-Doppler  radar  imaging,  in  [27].  Doppler  imaging 
of  a  planar  surface,  however,  does  not  seem  to  have  been  studied  in  its  own  right. 

The  purpose  of  this  paper  is  to  develop  the  theory  of  Doppler  imaging  for  the  case  of  a  single 
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sensor  moving  above  a  collection  of  scatterers  lying  on  a  surface.  In  section  2  we  set  forth  our  no¬ 
tation  and  a  model  for  the  received  radar  signal.  In  section  3  we  specialize  to  the  case  of  a  at  earth 
and  constant-velocity  ight  path.  For  this  case,  we  develop  an  approximate  image  reconstruction 
algorithm. 


2.2  Radar  data 


The  basic  physics  behind  radar  detection  is  simple  to  understand.  A  time-varying  voltage  sin(t )  £ 
C  is  passed  to  a  transmitting  antenna  where  it  activates  a  radiating  electromagnetic  eld  Em(x,t )  £ 
C3  de  ned  at  spatial  position  x  £  M3  and  time  t.  This  radiating  eld  obeys  the  vector  wave  equa¬ 
tion  and  interacts  with  radar  targets  by  inducing  time-varying  charge  distributions  j  (x,  t)  on  and 
within  their  support.  In  turn,  these  induced  currents  establish  a  response  electromagnetic  eld 
Esc  (x,  t )  which  radiates  to  a  radar  receiving  antenna  where  it  excites  an  echo  signal  voltage  ssc(t) 
that  is  fed  to  the  radar  receiver.  The  radar  system  then  uses  the  known  function  sm(t )  and  the 
measured  function  ssc(t )  to  estimate  various  properties  of  the  target  that  are  encoded  in  j  (x,  t )  (for 
example,  the  support  of  j  can  be  used  to  detennine  the  target’s  distance  and  bearing). 

The  details  of  the  physics  underlying  the  sequence  sm  — >  Ein  — >  j  — >  Esc  — >  ssc  is  usually 
quite  complicated  and,  in  practice,  various  simplifying  approximations  are  applied  to  keep  the 
analysis  tractable. 


2.2.1  The  scattered  field 


The  most  common  radar  echo  model  is  based  on  a  linearized  and  polarization-insensitive  scattering 
approximation  to  the  electric  or  magnetic  eld  integral  equations.  In  this  ad  hoc  approximation  to 
the  electromagnetic  scattering  problem,  we  consider  only  one  electric  eld  component,  namely  the 
one  for  which  the  antenna  is  designed  to  be  sensitive.  The  weak-scatterer  or  Born  approximation 
takes  the  induced  current  j  to  be  proportional  to  the  (scalar)  pEm,  where  p  is  the  target  reflectivity 
function,  and  Em(z ,  t  )  denotes  the  second  time  derivative  of  (one  component  of)  the  electric  eld 
incident  upon  the  target.  The  Born-approximated  scattered  eld  is  given  by 

E‘Bc(x,t)  ~-JJ  S(t  ~  ~  |*I/C0)PM  dt'dz  (2.1) 

where  Esfl{x,  t)  £  C  denotes  the  designated  component  of  the  scattered  electric  eld.  The  weak 
scatterer  approximation  is  well  documented  in  the  literature  (cf.,  [4,  5,  28,  39,  42,  45]  and  refer¬ 
ences  cited  therein)  and  will  not  be  further  motivated  here. 

We  assume  that  the  radar  emits  a  continuous  wave  at  the  x  ed  angular  frequency  u>0.  Then  the 
incident  eld  can  be  expressed  in  tenns  of  the  Green’s  function  for  the  three  dimensional  Helmholtz 
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equation: 


Ein{z,t)  =  EQ{Z^y)- 


a- iv0(t-\z-y\/c ) 


if  J  A  \  I  ’  V  *  / 

47r|z  —  y| 

where  y  denotes  the  phase  center  of  the  transmitting  antenna.  Substitution  into  equation  (2.1) 
yields 


e-i(o;o/c)(ct-|z-y|-|z-£c|) 

p(z)- — : - n - j - E0(z  -y)dz. 

\z  —  y\ \z  —  x\ 


Below  we  consider  only  the  case  where  the  transmit  and  receive  antennas  are  co-located,  have 
the  same  phase  centers  and  travel  along  the  ight  path  x  =  7  (t).  Moreover,  we  assume  that  the 
scattering  takes  place  in  a  thin  region  at  the  surface  z  =  C{zt),  where  zT  =  (z\ ,  z2).  Thus  we 
write  p(z)  =  p(zT)S(z  —  £(zr))-  We  also  introduce  the  notation:  R(zT,t )  =  C(zt)  —  7 (t); 
R(zT,t )  =  | R(zr,t)  | ;  and  R(zr,t)  =  R(zT,t)/R(zT,t).  The  eld  (2.3)  is  then 


EsBch(t),t)  =  J^r2  J  p(zTy 


2-iuo{t-2R(zT,t))/c) 

R2(zT,t ) 


E0(R(zT,t))dz7 


The  measured  signal  voltage  ssc(t )  arises  from  the  interaction  of  Esc  with  the  antenna  and  we 
can  write 


Ssc(t )  =  0 

'  1  J  (4tt)2 


'  e-iuj0(t-2R(zT,t))/c) 

p(zt )  - — r2(^Zt  ^ - W (R(zt,  t))  d zT  , 


where  W  is  a  weighting  factor  accounting  for  the  combined  transmit  and  receive  antenna  patterns 
[40]. 


2.2.2  Correlation  data 

The  radar  receiver  correlates  the  echo  signal  (2.5)  with  e~luJt,  which  is  a  frequency-shifted  version 
of  the  incident  signal  sin(t )  =  e~luJot,  over  a  nite  time  window;  in  other  words  it  measures  the 
windowed  Fourier  transform.  We  denote  by  i/)(t  —  r)  the  time  windowing  function  centered  at 
t  =  r;  we  write  the  (Fourier  transform)  frequency  as  u  =  ay//.  With  this  notation  the  radar  data 


rj(r,  p)=  ssc(t)eiulofJ,<'t~T'>'ijj{t  -  r)  d t 


^0 

(47r)2 


„—\u>o{t—2R(zr,t)/c) 

p(zt)—RE - - WWZT,  -  t )  dtdzT  . 

IX  yZjy ,  tj 


Then,  using  the  Taylor  expansion 


7(f)  =7(7)  +7  (r)(t-T)  + 
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to  write 


R(zT,t)  =  \C(zT)  =  ICOt)  -l(r)  -7(T)(t-T)H - | 

~  R{zt,t)  -R(zt,t)  -7(r)(t-r), 


in  the  exponent  of  (2.6),  we  obtain 


r/(r,  //) 


Wn  II  ,  .  e 


-itci0(t-2[ii(zT,T)-R(zT,T)-'y(r)(t-T)]/c) 


(4vr)2  JJ  p{Zt)~  K2(zt,t) 

x^(H(zr,t))e^‘-T)^(t  -  r)  dtd2T  • 


Finally,  in  (2.9)  we  make  the  change  of  variables  t'  —  t  —  r  to  obtain 

7 7(t,//)«  IJ  e-iwot'^+2^ZT’T>^T^  A(zT,t',T)dt'  p{zT)dzT 
where  t'  —  t  —  r  and 

p—iuJo{r—'2R[zT,T)/c)/jptf\ 

A '>  -  °  (W,r))a  ^  +  r))  • 


(2.8) 


(2.9) 


(2.10) 


(2.11) 


We  see  that  the  F  integration  results  in  an  approximate  delta  function  along  the  curve  formed 
by  the  intersection  of  the  earth  surface  with  the  constant-Doppler  cone  p  —  1  —  2 R(t,  Xt)  ■  7  (r)/ c. 
When  the  earth’s  surface  is  a  at  plane,  this  curve  is  a  hyperbola. 

We  use  the  dimensionless  parameter  (3  =  cuot'  to  write  A(zr.  (3,t ))  =  A(zT.  t/ .  r )  and 

<p{P,Zt,t,h)  =  P[l  ~p  +  2 R(zt,t)  -7 (t)/c]  (2.12) 

so  that  (2.10)  is 


»7(t,//)=  //  2l(zT,/3,r)e-W’z-^d/3p(zT)dzT. 


(2.13) 


2.3  The  Case  of  a  Straight  Flight  Path  and  Flat  Earth 


We  choose  coordinates  so  that  the  ight  path  is  along  the  z\  axis:  7(7")  =  (vr,0,H),  and  we 
assume  that  the  radar  is  operating  in  strip-map  mode,  meaning  that  the  antenna  beam  is  x  ed  and 
side-looking.  We  assume  that  the  antenna  beam  illuminates  the  region  \z\  —  vt\  <  e.  The  phase 
(2.12)  is  thus 


ip((3,ZT,T,p) 


p 


1  —  /i  + 


2v  Z\  —  VT 

c  \J (z\  -  VT )2  +  z\  +  H 2 


(2.14) 
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=  0 


1  —  /i  + 


2  v/c 


+  H 2 


(zi  -  vt)  +  0((z±  -  VT )2) 


where  in  the  second  line  of  (2.14)  we  have  done  a  Taylor  expansion  about  the  point  z\  =  vt.  We 
use  the  notation  p  =  p  —  1  and  use  0  for  the  approximate  form  of  the  phase: 


4>{P,zt,t,  //) 


2u/c  /zi 

V(^2/fT)2  +  1  Vlf 


(2.15) 


The  leading-order  contribution  to  rj  comes  from  the  critical  set  d(j)/df3  =  0,  which,  in  this  ap¬ 
proximation,  is  a  straight  line  in  the  p-r  plane.  If  instead  of  r  we  use  the  dimensionless  variable 
f  =  vt/H,  then  we  nd  that  the  straight  line  in  the  p-f  pane  has  f -intercept  Z\/H  and  slope 
—2 (v/c)((z2/H)2  +  l)-1/2.  We  note  that  zi  can  be  found  from  the  f-intercept,  and  z2  from  the 
slope. 


2.3.1  Image  formation 

We  form  an  image  by  backprojection: 

I(Vt)  =  Iff  eimyT-T^B(f3,yT,T,fi)df3ri(yT,jl)dlidT 


(2.16) 


where  B  is  to  be  determined  below. 


To  determine  the  degree  to  which  the  image  I  reproduces  p,  we  insert  (2.10)  into  (2.16)  and 
perform  a  stationary  phase  reduction  in  //  and  one  of  the  3  variables.  This  results  in 


I(yT)  =  /  K(yT,zT)p{zT)dzT, 


(2.17) 


where  the  point  spread  function  K  is  given  by 


2v 

K(yT,zT )  =  I/  exp  (  i/3— 


yi—vr 


Zl—VT 


VV22+H2  VZ2+H2 


A{zt,t) 


(2.18) 


x  B{/3,  yT ,  r,  y(yT,  t))i/}((3)  d(3dr 


and  where  p(yT,r )  =  —(2v/c)(yi  —  ur)(y2  +  H 2)  x/2.  We  would  like  to  choose  B  so  that 
K(yT,zT)  =  5(yT  -  zT)  =  (27t)'2/ 

The  leading-order  contributions  to  (2.18)  come  from  the  critical  points  of  (2.18),  which  are 
determined  by 


<90  y1-  vt 

0  =  oc 


Z\  —  VT 


d(3  ^/yf+JP  VW+H"2' 
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0 


1 


1 


(2.19) 


d(j) 

^  K  ^y(+JP~^f4  +  Jp' 

These  equations  have  solutions  only  when  (2/1, 2/2)  =  (zi,  dzz2),  which  implies  that  the  only  ar¬ 
tifacts  are  the  usual  left-right  ones.  (And  these  artifacts  can  be  avoided  by  a  judicious  choice  of 

W(R,t).) 


In  (2. 18),  we  expand  the  phase  using  the  identity 

f(zT )  -  f{yr )  =  (zT  -  Vt)  ■  [  V/(yT  +  Kzt  ~  yr))d\  «  (zT  -  Vt)  •  V/(yT)  (2.20) 


applied  to  /(yT)  =  (r/i  —  VT)(y?2  +  f/2)  x/2.  We  then  make  the  change  of  variables 


(P,  t)  -*•  £  =  (2 :/3u/c)V/(yr)  =  /?— 


1  -(2/1  -  i>r)  2/2  \ 

C  VV^  +  ^’  +  #2)3/2  ) 

where  the  superscript  “T”  denotes  transpose.  This  converts  (2.18)  into 


K{yT,zT)=  j  el(ZT  yTHA(zT,T)B(P,ii(yT,T),T,yT)ilj(P) 


d  (P,t) 


dt, 


(2.21) 


(2.22) 


where  now  (3  and  r  are  understood  to  refer  to  /3(£)  and  t(£).  The  Jacobian  in  (2.22)  is  the  reciprocal 
of 

1 

0 


<9(/3,r) 


4u2/3 


-(|/i  -  vr)y2 
{yl  +  ff2)3/2  (2/2  +  AT2)3/2 


vij2 


4v3/3y2 


c2(y2  +  H2)2  ' 


(2.23) 


From  (2.22),  we  see  that  £?  should  be  chosen  according  to 


B(P,  y,  t,  yT)  = 


x(P,  Ti  Vt) 


ae. 


9(0, t) 


(2vr  )2A(yT,T)^((3) 


(2.24) 


where  y  is  a  cutoff  function  that  prevents  division  by  zero  in  (2.24).  With  this  choice  of  B,  (2.22) 
becomes 


k(Vt,Zt)  =  A  f  fjMidS  (2.25) 

(2vr)2 

where  the  integration  is  over  the  set  Ey  of  £  swept  out  according  to  (2.21)  as  r  and  (3  range  over 
the  data  collection  region  for  the  point  yT. 


2.3.2  Resolution 

The  resolution  of  the  image  (2.16)  is  determined  by  the  region  Ey  of  integration  in  (2.25),  which 
is  the  set  in  (2.21)  as  r  and  (3  range  over  the  subset  of  the  data  collection  region  [rmin,  rmax]  x 
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[Anin,  Anax]  that  is  relevant  at  the  point  yT.  In  particular,  the  region  Ey  is  also  restricted  by  the 
beam  pattern  W :  the  maximum  X  of  x  =  \y\  —  vr\  is  the  distance  along  the  ight  path  for  which 
the  point  y  remains  in  the  beam.  Thus  x  is  in  the  interval  [— X,  A"]  where  X  =  vT  with  T  being 
half  the  persistence  interval. 

Since  (3  =  ix0t'  —  u>0(t  —  r),  the  length  of  the  interval  [Anin,  Anax]  is  2n  times  the  number  of 
cycles  in  our  time  window.  We  write  Anin  =  —  D  and  /3max  =  Q.  We  assume  that  the  antenna  beam 
is  directed  to  the  left  of  the  ight  path,  so  that  y2  >  0. 


The  boundary  of  the  region  5  is  formed  by  the  four  curves  Cq,  (An,  Cx,  C-x’- 

C±n  =  {(6,6 (x)):^  =  ±aQ,  ^  =  ±Qbx,  xe[-X,X}} 

C±x  =  {£(/3)  =  (a,±Xb)/3:l3e[-n,n]}  (2.26) 

where  we  have  written  a  =  2 (v/c)(y2  +  if"2)-1/2  and  b  =  —  2(v  /  c)y-2(y2  +  if2)-3/2. 

We  see  that  for  a  given  target  location  (y\ ,  y2),  the  curves  C±q  are  vertical  lines  (in  which  only 
6  varies).  Similarly,  the  curves  C±x  are  radial  lines  through  the  origin.  The  region  Ey  is  thus 
a  bowtie.  The  size  of  Ey  depends  not  only  on  the  target  location  (2/1, 2/2)  but  also  on  the  system 
parameters  v,  H ,  lu0. 


To  obtain  the  point  spread  function  (ambiguity  function)  for  a  particular  image  point  (0, 2/2)5 
we  calculate  the  right  side  of  (2.25).  We  write  p  =  zT  —  yT  in  (2.25)  and  change  variables  from 

|  =  (a,  bx)(3  to  (x,  (3),  obtaining 


Ky(P) 


1 


eX 


3 i(a,bx)-p0 


l-X 


9£ 

9(/3,  x) 


dx  d (3. 


(2.27) 


The  Jacobian  \d£/d(f3,x)\  is  easily  found  to  be  equal  to  \ab/3\;  K  can  be  calculated  as  follows. 


Ky(p)  = 


-1 


(2vr)2 

pi (a,bX)-pQ  _ 


cX 


-a  J -x 


j(a,bx)  ppabj3  da;  dp  _ 


-O  rx 


j(a.hx)-pf) abp  da;  d/3 


'0 


i-x 


pi (a,-bX)-pQ  _  p-~i(a,bX)-pfl  ^  p— i(a,— b.Y)-pf2 

+ 


(27r)2ip2  V  i (a,bX)-p  i(a,  —  bX)  ■  p  i (a,bX)  ■  p  i (a,-bX)-p 

a  f  sm2((a,bX)  ■  ptl/2)  sin2((a,  —  bX)  ■  p  0/2) 


n2p2 


(a,  bX)  p 


(a,  —bX)  p 


.(2.28) 


Down-range  resolution. 

Ky(0,P2)  =  ( 

1T2P2  \ 

where  we  have  used  sin2  A 
resolution  is 


Down-range  resolution  is  obtained  by  setting  pi  =  0  in  (2.28): 

sin2(06Xp2/2)  sin2(— QbXp2/2)  \  sin2(D6Xp2/2) 

bXp2  —bXp2  j  *  (VlbXp2/2)2 

=  sin2 (—A).  From  the  right  side  of  (2.29),  we  see  that  the  down-range 


4t r  _  2vr c(y2  +  H2f'2 
QbX  ~  nxvy2 


(2.30) 


13 


where  the  resolution  is  de  ned  as  the  width  of  the  central  lobe  of  (2.29).  We  see  that  the  down- 
range  resolution  is  improved  by  taking  a  larger  synthetic  aperture  2X;  this  is  consistent  with  our 
expectation  that  a  single  Doppler  measurement  will  provide  no  range  information.  Range  informa¬ 
tion  is  obtained  only  from  combining  measurements  across  the  synthetic  aperture. 

For  uq  =  2n  x  100  MHz,  v  —  8  m/s,  a  0.5  second  time  window,  and  y2  =  400,  the  down-range 
resolution  is  roughly  A p2  ~  80  m. 


Cross-range  resolution.  Cross-range  or  azimuthal  resolution  is  obtained  by  calculating  the  p2 
0  limit  of  (2.28) ;  this  can  be  done  by  means  of  l’Hospital’s  rule: 


lim  Ky(p1,p2)  oc 

P2  *0 


sin(f2opi/2)  cos(f2api/2)  sin2(f2api/2) 


Qapi/2 

f lap,  'Sill(n“Pl) 
1  sin(Qapi) 

Qapi  ttapi 
1  sin  (Dap  i) 
ttapi  Qap\ 


2(QaPl/2y 


siir(Hapi/2) 


Qapi 


flapi/2 

1  —  cos(Qapi) 


sin(Oapi) 
[flapi  —  tan(Oapi/2)]  . 


(2.31) 


The  width  of  the  central  lobe  of  Ky{p\ ,  0)  is  an  estimate  of  the  cross-range  resolution  and  the 
rst  zero,  which  arises  from  the  tenn  in  brackets  in  the  last  line  of  (2.31),  involves  solving  the 
transcendental  equation  a  —  tan(a/2)  =  0,a^0.  We  obtain  a  =  2.33  . . .,  and  so  the  cross-range 
resolution  is 


A  pi  ~ 


2  x  2.33 
f la 


2.33  c(y\  +  H2)1'2 
Vtv 


(2.32) 


For  the  same  parameter  values  given  above,  the  cross-range  resolution  is  roughly  Api  «  100 
m. 


We  see  that  both  down-range  and  cross-range  resolution  is  improved  by  having  a  larger  number 
of  wave  cycles  within  the  time  window  (i.e.,  by  increasing  the  frequency  cj0  or  using  a  longer  time 
window)  and  increasing  the  ight  velocity  v.  Resolution  is  worse  for  points  more  distant  from  the 
ight  path. 


2.3.3  Numerical  Implementation 

A  full  implementation  of  the  scheme  in  section  2.3.1  is,  computationally,  very  expensive.  Equation 
(2.16)  is  a  triple  integral  with  kernel  B  given  by  equation  (2.24).  Moreover,  evaluation  of  this 
ltering  kernel  requires  calculation  of  A  (equation  (2.11))  and  the  Jacobian  (equation  (2.23)).  Even 
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with  the  use  of  look-up  tables,  the  computation  requirements  associated  with  such  spatially-varying 
kernels  can  be  daunting. 


Instead,  we  make  use  of  our  previous  observation:  the  leading-order  contribution  to  rj  comes 
from  the  critical  set  d<f>/df3  =  0  which,  in  this  approximation,  is  a  straight  line  in  the  /i-r  plane. 
Consequently,  an  alternative  approach  can  be  based  on  identifying  the  straight  line  components  in 


An  ef  cient  line-detection  algorithm  makes  use  of  the  Radon-Hough  transfonn  that  integrates 
over  all  possible  lines: 


H{rj}(d,d) 


rj(r,  jl)8(d  —  t  cos  9  —  /isinf?)  drd/i. 


(2.33) 


The  transform-image  of  77  yields  a  detected  (matched)  line’s  perpendicular  offset  from  the  origin  d 
and  the  angle  9  of  this  offset.  The  offset  d  and  angle  6  are  related  to  the  slope  m  and  r-intercept 
b  via  d  =  b  cos  6  and  m  =  tanff  Because  r  has  the  dimensions  of  time,  m  has  the  dimensions 
of  reciprocal  time.  These  dimensional  quantities  are  related  to  the  dimensionless  r-intercept  b  and 
slope  m  of  the  line  in  the  /i-f  plane  via  b  =  vb/ H  and  fh  =  Hm/v.  According  to  our  analysis 
above,  we  can  therefore  form  the  estimates 

Zi  =  — —  and  cote=—  (=^)  J (f )  + 1  (2.34) 

cos  9  v  v 

_  _  1  /  4-m3  cot2  9  7 

^  ^2  ~  h  V  “^2  1 


to  locate  the  scattering  center.  Our  expedient  imaging  scheme  is  therefore  performed  as  a  two- 
step  process:  1)  form  the  Radon-Hough  transform  H{r]}(d,  9)  of  the  correlation  data;  and  2)  map 
(resample)  according  to  equation  (2.34). 


Figure  3.7  shows  the  (magnitude  of  the)  Doppler  frequency  response  of  three  isolated  equal- 
strength  scatterers  as  a  function  of  time.  (The  difference  in  relative  strength  of  each  curve  is  due  to 
a  difference  in  scatterer  range.)  The  simulated  data  were  constructed  under  the  model  assumptions 
of  a  straight  ight  path  for  the  radar  with  transmit  frequency  of  cu0  =  100  MHz  and  relative  velocity 
of  8  m/s.  Each  data  point  was  determined  from  a  0.5  sec  modeled  measurement  interval  (according 
to  equation  (2.13)). 


Figure  3.8  was  formed  from  the  data  of  Figure  3.7  by  resampling  the  Radon-Hough  transform 
according  to  equation  (2.34).  The  gure  displays  the  magnitude  of  the  simple  transfonn  and  no 
attempt  has  been  made  to  conect  for  range-dependent  magnitude  enors — although  this  next  image 
processing  step  could  be  easily  made  once  the  individual  scatterer  locations  have  been  determined. 
Figure  3.8  could  be  improved  by  image-processing  techniques  such  as  high-pass  ltering  or  thresh¬ 
olding.  We  leave  this  topic  for  the  future. 


We  note  that  Figure  3.8  shows  a  superposition  of  three  point-spread  functions  for  our  expedient 
imaging  scheme. 
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Figure  2.3:  The  Doppler  return  from  three  point  scatterers. 


-1000  -500  0  500  1000 

Zi  (meters) 

Figure  2.4:  The  (resampled)  Radon-Hough  transform  of  the  data  in  Figure  3.7.  The  true  locations 
are  marked  by  o;  the  ight  path  is  along  the  z\  axis. 

2.4  Conclusions  and  Future  Work 


We  have  developed  a  mathematical  model  for  the  radar  return  signal  resulting  from  a  transmitted 
CW  waveform  in  the  case  of  an  antenna  moving  above  a  scattering  surface.  For  the  case  of  a  ight 
path  with  constant  velocity  over  a  at  surface,  we  have  developed  two  imaging  algorithms.  Finally, 
we  have  displayed  the  results  of  numerical  simulation  for  a  simple  case. 

An  important  implementation  issue  that  we  did  not  address  concerns  the  dif  culty  of  collecting 
these  data  from  monostatic  systems  (in  which  the  transmitter  and  receiver  are  co-located).  Since 
the  incident  wave  is  constantly  transmitting  at  frequency  cuo,  any  echo  wave  with  zero  Doppler 
shift  will  be  undetectable  in  the  correlation  receiver — the  echo  signal  strength  will  be  too  small  to 
compete  with  that  of  the  transmitter.  It  is  conceivable  that  a  judicious  choice  of  antenna  weighting 
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function  W (R)  combined  with  frequency  domain  ltering  could  be  used  to  “notch-out”  this  region 
of  the  data.  Such  a  notch  would  need  to  be  suf  ciently  narrow  that  only  a  small  fraction  of  the 
available  data  are  excluded  and,  for  low  frequency  systems,  the  required  narrowness  of  this  notch 
might  be  impracticable.  Of  course,  this  problem  is  not  generally  manifest  in  bistatic  systems  for 
which  the  transmitter  and  receiver  are  widely  separated  and  our  analysis  can  be  readily  modi  ed 
to  include  such  radar  con  gurations. 

Much  more  work  remains  to  be  done,  especially  for  cases  in  which  the  geometry  is  more 
complicated.  In  addition  we  leave  for  the  future  the  problem  of  developing  fast  numerical  imple¬ 
mentations  [32]  of  the  general  imaging  fonnula. 
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Chapter  3 


Imaging  that  exploits  mulitpath  scattering 

This  work  is  joint  with  Bob  Bonneau;  the  theory  has  been  published  in  [9]. 


3.1  Introduction 


This  work  is  motivated  by  the  problem  of  radar  imaging  in  urban  and  forested  settings.  Trees 
contain  branches,  which  could  cause  radar  re  ections;  thus  an  object  on  the  ground  could  be  il¬ 
luminated  not  only  by  waves  traveling  a  direct  path  from  the  antenna  to  the  object,  but  also  by 
waves  that  have  been  re  ected  from  branches  onto  the  object.  Similarly,  imaging  in  urban  settings 
involves  multipath  scattering. 

Multipath  scattering  confuses  standard  imaging  methods.  However,  multiply-re  ecting  waves 
illuminate  the  object  from  a  variety  of  directions;  if  these  waves  can  be  accounted  for  and  used  in 
the  imaging  process,  they  might  improve  the  image  of  the  object.  This  paper  provides  an  imaging 
method  that  exploits  such  multipath  scattering.  We  nd  that  artifacts  can  arise  in  certain  situations; 
we  show  how  to  avoid  them  and  analyze  the  improvement  in  resolution  due  to  the  multiply  scattered 
waves.  We  carry  out  the  details  of  the  analysis  explicitly  for  the  case  of  a  single  extra  scattering 
path;  the  extension  to  more  paths  is  expected  to  follow  along  similar  lines. 

We  consider  the  case  in  which  the  multipath  scattering  is  due  to  the  presence  of  scattering 
“centers”  in  the  foreground.  These  scattering  centers  we  model  as  point  scatterers.  Such  point 
scatterers  are  commonly  used  to  model  scattering  from  comers  and  edges.  Using  point  scatter¬ 
ers  has  the  advantage  that  an  exact  closed-form  solution  is  available  [37]  for  multiple  scattering 
between  such  scatterers. 

For  scattering  from  the  target  or  scene  of  interest,  we  use  the  Born  (single-scattering)  approx¬ 
imation.  This  approximation  neglects  multiple  scattering  within  the  target  and  multiple  scattering 
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Figure  3.1:  Scattering  geometry 


between  the  target  and  its  environment;  this  approximation  is  ubiquitous  in  radar  work. 

In  the  paper  we  discuss  explicitly  the  radar  case,  but  the  analysis  applies  equally  well  to  sonar 
and  ultrasound  imaging. 

The  paper  is  organized  as  follows.  In  section  3.2  we  develop  a  mathematical  model  for  the 
measured  signal.  We  show  that  this  signal  is  of  the  form  of  a  Fourier  integral  operator  applied  to 
the  scene.  In  section  3.3  we  give  a  method  for  producing  an  image,  and  show  that  the  image  has 
the  desired  properties.  In  section  3.4  we  discuss  the  resolution  obtained  from  this  method. 


3.2  The  Mathematical  model 


The  mathematical  model  involves  a  number  of  ingredients:  1)  a  model  for  wave  propagation  in 
free  space,  2)  a  model  for  the  antenna,  3)  a  model  for  multiple  scattering  from  point  scatterers  in 
the  foreground,  and  4)  a  (linearized)  model  for  scattering  from  the  scene, 


3.2.1  A  model  for  the  wave  propagation 

Each  component  of  the  electromagnetic  eld  propagates  in  free  space  according  to  the  scalar  wave 
equation: 

V2w  —  c~2ii  =  0, 

where  the  dots  denote  differentiation  with  respect  to  t. 


(3.1) 


Equation  (3.1)  is  not  a  perfect  model  for  electromagnetic  scattering:  when  electromagnetic 
waves  encounter  materials,  there  is  coupling  between  the  different  vector  components  of  the  elec¬ 
tromagnetic  eld.  By  using  (3.1),  we  are  ignoring  such  coupling,  and  we  are  thus  ignoring  polar¬ 
ization  effects. 

We  will  use  capital  letters  for  frequency-domain  quantities,  which  are  related  to  time-domain 
quantities  by  the  Fourier  transform 

U{uj,x)  —  —  f  elu>tu(t,x)dt.  (3.2) 

2vr  J 

We  write  k  =  u/c0. 

We  consider  the  case  in  which  N  scattering  “centers”  are  present  in  the  foreground  of  the  scene 
we  wish  to  image.  These  scattering  centers  we  model  as  point  scatterers  5zj(x )  =  S(x  —  z] ) ,  j  = 
1,2, ...  ,N.  If  we  assume  that  the  eld  is  due  to  an  isotropic  point  source  at  position  y  and 
time  t  —  0,  the  corresponding  differential  equation  is 

N 

v lgN(t,x,y )  -  Co2gN(t,x,y )  -  ^ yiSzi(x)gN(t,x,y)  =  -Sy(x)S(t)  (3.3) 

i=i 

Here  the  y’s  are  the  strengths  of  the  point  scatterers.  An  explicit  expression  for  gjy  is  given  in  the 
next  section. 

Behind  the  cloud  of  scattering  centers  is  the  scene  or  target  we  wish  to  image.  We  assume  that 
scattering  in  the  scene  of  interest  is  due  to  a  perturbation  of  the  wave  speed,  which  we  denote  by  q. 

Whether  we  can  form  a  three-dimensional  image  of  a  volume  or  merely  a  two-dimensional 
image  of  a  surface  is  detennined  by  the  number  of  degrees  of  freedom  in  our  measured  data.  In 
the  case  of  a  (  x  ed)  steerable  array  antenna,  for  example,  we  measure  a  time-varying  eld  for  each 
steered  direction.  If  the  beam  can  be  steered  in  two  directions,  we  have  three  degrees  of  freedom 
in  the  data  and  can  expect  to  make  a  three-dimensional  reconstruction. 

In  the  strip-map  synthetic-aperture  radar  case,  however,  we  measure  a  time-varying  signal  from 
each  point  along  a  one-dimensional  curve;  the  data  then  contains  two  degrees  of  freedom  and  we 
expect  to  make  only  a  two-dimensional  image  of  objects  sitting  on  a  known  surface.  In  this  case, 
we  denote  the  known  surface  by  {a?  =  ij){xT)  :  xT  e  M2}),  where  xT  =  (aq,  x2),  and  we  write 
the  wave  speed  perturbation  as  q(xT)S(x  —  ij)(xT)). 

In  what  follows  we  consider  for  simplicity  the  three-dimensional  case;  for  the  two-dimensional 
case  we  would  simply  replace  x  by  xT  and  q(x)  by  q(xT)S^(x  —  i/>(c cT)). 

The  equation  we  consider  can  thus  be  written 

N 

V2w  —  Cq  2u  —  Hi5ziU  —  Cq  2gii  =  0  (3.4) 

i= 1 

In  (3.4),  q  is  the  quantity  we  wish  to  reconstruct. 
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3.2.2  A  model  for  the  antenna 


We  model  the  antenna  as  an  array  antenna  in  which  each  element  is  an  isotropic  radiator.  At 
frequency  uj,  the  eld  at  x  emanating  from  the  element  at  y  is  G0(cu,  \x  —  y\),  where  G0(cj,r)  = 
(47rr)‘  1  exp(i kr).  To  model  the  eld  u  from  an  extended  antenna,  centered  at  y,  over  which  the 
derivative  of  the  current  density  is  j(t,y'),  where  y’  denotes  a  vector  from  the  center  y  to  an 
arbitrary  point  on  the  antenna,  we  simply  integrate: 

e~luJt  f  G0(u,\x  —  (y  +  y')\)J(u,y')dy'dw.  (3.5) 

J  antenna 

(Here  J  denotes  the  Fourier  transform  of  j.) 

A  steered  antenna  beam  can  be  modeled  by  including  phasing  in  j.  For  example,  for  a  beam 
steered  in  direction  y,  one  can  choose  =  exp(-ic jy  ■  y’)\ 'antenna (j/').  where  Xantenna 

denotes  a  function  that  is  one  on  the  antenna  and  zero  off  it. 

In  what  follows,  we  consider  only  the  eld  Go  emanating  from  a  single  antenna  element,  with 
the  understanding  that  these  elds  can  be  assembled  to  model  any  desired  antenna. 

Reception  by  the  antenna  is  modeled  by  integrating  the  scattered  eld  over  the  antenna,  again 
perhaps  with  some  weight  for  steering.  For  simplicity  we  consider  only  the 


3.2.3  Multiple  scattering  from  point  scatterers 

For  a  time-harmonic  incident  wave  Um,  the  frequency-domain  eld  Usc  scattered  from  N  “point” 
scatterers  can  be  obtained  from  the  Foldy-Lax  [45]  equations  together  with  the  assumption  that 
the  scattered  eld  from  an  individual  “point”  scatterer  is  proportional  to  the  free-space  Green’s 
function  G0  [37]: 

N 

Usc(uu,x)  =  Gq(u,  \x  —  zj\)fijUj(u! ,  zj)  (3.6) 

j= i 

Uj(u,x)  =  Uin(u,  x)  +  Gq(u,  \x  -  z'DiijUjiu,  zl),  j  =  1, 2, . . . ,  N, 

*7 6? 

(3.7) 

Equation  (3.6)  says  that  the  scattered  eld  is  the  sum  of  the  elds  scattered  from  each  scatterer; 
moreover,  the  eld  scattered  from  the  jth  scatterer  is  proportional  to  the  eld  Uj  that  is  incident 
upon  the  jth  scatterer.  Equations  (3.7)  say  that  the  yth  local  incident  eld  is  the  overall  incident 
eld  plus  the  eld  scattered  from  all  the  other  scatterers.  If  the  scattering  strengths  yi,  y2,  ■  ■  ■ ,  I^n 
and  positions  z1,  z2, . . . ,  zN  are  known,  the  equations  (3.7)  can  be  solved  for  the  Uj;  then  the  total 
eld  U  =  Umc  +  Usc  can  be  found  from  (3.6). 
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The  total  eld  U  satis  es  the  “background”  differential  equation 

N 

V2U  +  k2U  +  J2  =  0  (3.8) 

i— 1 

where  k  =  u/cq.  We  note  that  the  sense  in  which  a  eld  of  the  form  (3.6)  satis  es  (3.8)  requires 
an  extension  of  the  traditional  distributional  de  nition  of  the  delta  function;  its  domain  must  be 
extended  to  include  functions  with  1/r  singularities.  A  thorough  discussion  of  this  issue  can  be 
found  in  [2], 

We  will  write  Gn  for  the  full  eld  U  in  the  case  when  Umc  =  G0. 


Example:  A  single  point  scatterer 

For  a  single  point  scatterer  located  at  position  z,  the  scattered  eld  is  simply 

Um(x)  =  k2G0{\x  -  z\)/iUm(z)  (3.9) 

The  corresponding  time-domain  scattered  eld  is 

um(t,x)—  [  e~lwtG0(oj,  \x  —  z\)yUin(uj,  z)du>.  (3.10) 


We  obtain  the  one-point-scatterer  “background”  Green’s  function  G\  by  taking  Um 
Gi(u,x,y')  =  G0(cu,  \x  -  y  |)  +  G0(u,  \x  -  z\)yG0(uj,  \z  -  y'\) 

We  will  denote  the  second  term  on  the  right  side  of  (3. 1 1)  by  Gf . 


G  o  : 

(3.11) 


Example:  A  pair  of  point  scatterers 


In  the  case  of  two  “point”  scatterers,  equations  (3.7)  are 


EM*)  =Uin(x)  +  G0(\x-z1\)y2U2(z2) 
EM*)  =  Uin(x)  +  G0(\x  —  z‘2\)yiUi(z1) 


Evaluating  (3.12)  at  z1  and  (3.13)  at  z~  gives  rise  to  the  system  of  equations 

(  1  -^Go(L)  \  (  Ui(zl)  \_(  Uinc{zl)  \ 

V  ~t*iG0(L)  1  )  \  U2(z>)  J  ~  {  Uinc(z2)  J 

where  L  =  \z2  —  z1].  These  equations  have  the  solutions 


UjW 


- 


Uinc(zi )  +  yfG0{L)Uinc{zj') 

1  —  /J,i112Gq(L) 


3  =  1,2, 


(3.12) 

(3.13) 


(3.14) 


(3.15) 
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where  j'  —  2  if  j  —  1  and  j'  =  1  if  j  =  2.  Using  (3.15)  in  (3.6)  yields 


Us' 


J  =  1 


33 


iv. 


(/inc(^)+<UyG0(L)(/il 
1  —  hi/12Gq(L) 


(3.16) 


Equation  (3.16)  has  a  clear  physical  interpretation  if  we  consider  the  denominator  to  be  the 
sum  of  a  geometric  series: 


2  OO 

UK(x)  =  X)  [G°(l*  -  *'l) 

j=l  n=0 

oo 

+  G0(|*  -  ^'|)^Go(L)  ^[/^GgtL)]"^;^ 

n=0 


(3.17) 


The  n  —  0  term  on  the  rst  line  of  (3.17)  corresponds  to  the  incident  wave  scattering  once  from 
z3 ;  the  n  —  0  tenns  on  the  second  line  corresponds  to  the  incident  wave  scattering  once  from 
z3  and  once  from  z] .  The  n  —  1  tenn  on  the  rst  line  corresponds  to  initial  scattering  from  z3 , 
then  from  zJ  ,  and  then  from  zJ  again.  The  n  —  1  term  on  the  second  line  corresponds  to  initial 
scattering  from  z3  and  two  bounces  off  z3 .  The  terms  corresponding  to  larger  values  of  n  have 
similar  interpretations. 


In  any  physical  problem,  some  energy  loss  occurs  with  each  bounce  (modeled  by  the  /i’s  being 
less  than  one),  so  that  only  a  few  tenns  in  the  series  are  relevant. 


The  two-point-scatterer  “background”  Green’s  function  G2  is  found  by  taking  Umc  =  G0: 


2 

G2(w,x,y')  =  G0(u,  33,  y')  +  ^  G0(u,  \x  —  z3\)nj 

3= 1 

Go(uj,  | z3  —  y'\)  +  Hj'Go(ui,  L)Gq(u) ,  | z3  —  y'\) 

1  —  fii^G^u),  L ) 


(3.18) 


N  point  scatterers 


When  N  scatterers  are  present,  (3.7)  is  a  system  of  N  equations  that  can  be  solved  by  Cramer’s 
rule,  which  results  in  a  complicated  but  closed-form  expression  for  the  solution.  This  expression 
has  a  denominator  containing  the  determinant  of  coef  cients;  this  determinant  has  an  expansion 
that  allows  for  a  multipath  interpretation  similar  to  the  one  above.  As  before,  only  a  limited  number 
of  terms  need  to  be  retained. 
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3.2.4  The  model  for  scattering  from  the  target 


The  total  eld  G(u>,y,y')  at  y  due  to  the  antenna  element  at  y'  is  equal  to  the  sum  of  the  following 
elds:  a)  the  free-space  eld  G0(tu,  \y  —  y'\)  emanating  from  the  antenna  element,  b)  the  eld 
Gn  scattered  from  the  N  point  scatterers  in  the  foreground,  and  c)  the  eld  Gsc  scattered  from  the 
target  q(x). 

For  Gsc,  we  use  the  Born  approximation  or  single-scattering  approximation  to  model  the  scat¬ 
tered  eld. 

The  Bom  approximation  in  this  case  is 

y ,  y)  =  -  I  Gn(uj,  y,  x)q(x)GN(uj,  x,  y')u2dx.  (3.19) 

The  Born  approximation  makes  the  mapping  from  q  to  usc  linear,  but  it  is  not  necessarily  a  good 
approximation.  Another  linearizing  approximation  that  can  be  used  for  re  ection  from  smooth 
surfaces  is  the  Kirchhoff  approximation,  in  which  the  scattered  eld  is  replaced  by  its  geometrical 
optics  approximation  [8]  [25].  Here,  however,  we  consider  only  the  Bom  approximation. 

The  corresponding  time-domain  eld  is 

<7b(G  V:  y')  —  J  e~luJtGN(u,  y,  x)q(x)GN(uj,  x,  y')uj2dxdw.  (3.20) 

We  note  that  this  eld  is  of  the  form 

sfB{t,y,v')=  (3-21) 

j/G  {paths} 


where 

Fj[q](t,y,y')  =  j  e~1^[t~TjM’x)]aj(u,y,y',x)dwq(x)dx  (3.22) 

where  t3  denotes  the  travel  time  along  path  j  and  where  aj  contains  the  geometrical  spreading 
factors,  multiples  of  47T,  and  (the  Fourier  transform  of)  the  waveform  sent  to  the  antenna. 

The  total  eld  is  G  =  G0  +  GN  +  Gsc  «  G0  +  GN  +  Gs£. 


3.3  Image  formation 


We  outline  rst  the  strategy  for  the  general  case  of  N  point  scatterers,  then  carry  out  the  detailed 
analysis  for  the  case  of  a  single  point  scatterer. 
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3.3.1  General  strategy 


We  assume  the  the  source  location  y'  and  receiver  location  y  are  known;  thus  G0(ou,  y.  y')  can  be 
subtracted  from  the  received  eld.  This  leaves  Gr  =  G  —  G0,  which  we  consider  to  be  the  data  in 
the  image  formation  process. 

The  next  step  in  the  image  formation  process  is  to  identify  the  foreground  scatterers.  This  can 
be  done  from  the  early-time  part  ( Gn )  of  the  signal,  because  the  scatterers  are  assumed  to  be  in 
front  of  the  scene  of  interest.  Identi  cation  of  the  locations  z]  and  strengths  /j7  from  Gn  can  be 
done  in  a  number  of  ways.  One  approach  is  to  use  optimization,  in  which  one  nds  the  z’s  and  //’s 
that  minimize 

___ ;  __  ||  /^measured  /^calculated  1 1  / o  ^  o  \ 

mm  W^n  ~  (-ttv  II-  1  j.zj) 

zJ’H 


Another  approach  for  nding  the  locations  is  to  use  Devaney’s  MUSIC  algorithm  [17]. 
Although  the  treatment  in  [17]  is  based  on  the  Bom  approximation,  in  fact  Devaney’s  approach  ap¬ 
plies  also  to  the  multiple-scattering  case:  determining  the  locations  of  the  point  scatterers  depends 
only  on  the  fact  that  (3.6)  is  a  linear  combination  of  the  functions  Go(u,  \x  —  zj\). 

Once  the  locations  and  strengths  of  the  foreground  scatterers  are  known,  then  Gn  is  known  and 
can  be  subtracted  out.  This  leaves  Gsc,  from  which  we  form  an  image  by  ltered  backprojection: 

I(p)  -  ~  Y1  J  e'Ut-rj(y’I,,’pA(w>  p,  y,  y,  y!)dtdydy'  (3.24) 

j£  {paths} 

where  the  Iters  bj  are  determined  below.  We  note  that  in  (3.24),  the  paths  are  known  because  the 
foreground  scatterers  are  known.  We  see  below  that  we  must  take  precautions  to  avoid  artifacts  in 
the  image. 

In  (3.24)  the  integration  over  y  and  y'  indicates  that  we  sum  over  all  the  data.  As  written, 
(3.24)  is  appropriate  for  a  collection  of  isotropically-radiating  sources  and  receivers;  in  the  case 
of  a  single  steered  antenna,  we  would  integrate  instead  over  the  steering  vectors  [//  in  the  notation 
of  (3.5)]  for  the  transmitted  and  received  elds.  We  see  below  that  in  order  to  avoid  artifacts,  we 
backproject  only  along  paths  that  include  a  direct  path  to  (or  from)  the  scatterer. 

We  illustrate  the  imaging  process  and  its  analysis  for  the  case  of  one  foreground  point  scatterer. 


3.3.2  Case  of  a  single  point  scatterer 

For  the  case  of  a  single  point  scatterer  at  position  z,  the  Bom-approximated  eld  Gb  is  of  the  form 

9B(t,y,y')  =  gi(t,y,y')  +  [  e^tG1(u,y,x)q(x)G1(u,x,y')uj2dudx 
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(3.25) 


=  9o  +  9?  + J  e"Swt(G0  +  Glc)q(G0  +  Gf)u2dudx 
—  9o  +  9?  +  (^i  +  F2  +  F3  +  Fa)  [q] 

where  G\  is  given  by  (3.11),  where  Gf  =  G i  —  G0,  and  where  the  operators  f)  are 


Fi[q](t,y,y') 
F2[q}{t,y,y') 
F3[q}{t,y,y') 
F-Mlit,  y,  y') 


J  e~luJtG0(u,  y,  x)G0(u,  x,  y')q{x)dx 
J  e~luJtGf(uj,  y,  x)G0(u>,  x,  y')q(x)dx 
J  e~lujtG0(u y,  x)Gs^(u,  x,  y')q{x)dx 
J  e_lwtGf  (w,  y ,  x)G\c(lj,  x,  y')q(x)dx 


(3.26) 


The  different  F’s  correspond  to  different  scattering  paths:  F  corresponds  to  the  direct  path  from 
y'  to  the  target  to  y:  F2  corresponds  to  the  path  for  which  a  wave  leaves  y' ,  scatters  from  the 
foreground  scatterer  at  z  onto  the  target,  and  takes  a  direct  path  back  to  r/;  etc.  The  corresponding 
travel  times  r:)  for  these  paths  are 


Ti(y,y\p)  = 

\y  -  x  +  \x  -  y'\ 

F  {y,y,p)  = 

\y  —  z  +  \z  —  x  +  \x  —  y'\ 

T3{y,y',p)  = 

\y—x\  +  \x  —  z\  +  \z—-  y'\ 

n(y,y',p)  = 

\y  -  z  +  2\x  —  z\  +  \z  —  y'\ 

(3.27) 

Identification  of  scatterer  in  the  foreground 


We  assume  that  the  foreground  point  scatterers  are  closer  to  the  sensor  than  the  object  q  and  that 
therefore  the  early-time  part  of  (3.25)  consists  only  of  the  term  gi(t,  y,  y').  Since  we  know  y  and 
y',  we  can  subtract  out  g0  from  y, ,  leaving 


9si(t,  y,  y') 


J  e  10JtG0(u,\y  -  z\)yG0(u,\z  -  y'\)du 

r  e-iw[t-|y-z|-|z-y'|] 

^  J  (4vr)2|r/  -  z\\z  -  y' 


(3.28) 


In  the  case  of  a  single  point  scatterer,  there  can  be  no  multiple  scattering,  which  implies  that  the 
eld  g\  is  the  same  as  its  Bom  approximation.  We  can  thus  form  an  image  of  the  scatterer  2  by 
backprojection  as  described  below  for  the  i  =  1  case. 


Backprojection 


We  form  the  image  /  by  means  of  (3.24).  In  the  analysis  below,  we  replace  gsc  by  g^.  If,  in  (3.24), 
we  naively  backproject  along  all  possible  paths,  we  will  see  that  some  paths  cause  artifacts  in  the 
image. 
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Using  gsc  «  ^  Fj  [q]  [from  (3.22)]  in  (3.24)  results  in  an  equation  of  the  fonn 

Up)  = 

i=l  j= 1 

where  the  kernel  K  is  the  imaging  point-spread  function.  If  we  had  K(p,  x )  =  S  ip  —  x ),  then 
the  image  I  would  be  perfect;  we  want  to  determine  the  bj  so  that  K  comes  as  close  as  possible  to 
being  a  delta  function. 

The  contribution  to  K  from  BiFj  is 

Kid(p,x)  =  f  jUt-rM^b^p^y^y') 

xe_la/ d-TdyF 'xiaj(co' .  y,  y  ,  x)dujdu)' dtdydy'  (3.30) 

In  (3.30)  we  carry  out  the  t  and  to'  integrations,  obtaining 

Kid(p,x)  =  2vr  J  e™(Tilv'v''*)~Ti{y'y''p))hi(u,p,y,y,)a ,j(uj,y,y' ,x)dujdydy'  (3.31) 

In  order  for  the  K  to  be  a  close  approximation  to  a  delta  function,  we  would  like  the  diagonal  tenns 
Kl  t  themselves  to  be  good  approximations  to  delta  functions,  and  we  would  like  the  off-diagonal 
terms  Kij,  i  f  jj  to  be  zero  or  to  contribute  only  higher-order  tenns.  To  determine  whether  this  is 
the  case,  we  analyze  each  term. 

In  each  case,  the  analysis  is  similar:  we  use  the  method  of  stationary  phase  (see  appendix)  to 
determine  the  leading-order  contributions.  Analysis  of  the  critical  points  of  the  phase  detennines 
the  locus  of  points  p  that  will  appear  in  the  image  due  to  a  scatterer  located  at  x.  We  would  like 
the  the  critical  conditions  to  imply  that  p  =  x  ;  if  this  is  not  the  case,  the  other  possible  solutions  p 
tell  us  what  artifacts  will  appear  in  the  image  due  to  a  scatterer  at  x. 

Analysis  of  the  critical  conditions  detennine  where  in  the  image  a  scatterer  at  x  is  positioned; 
the  coef  cients  i),  detennine  the  amplitude  of  the  image  at  x.  The  coef  cients  1),  we  determine 
from  the  diagonal  terms  Kt  l.  In  these  diagonal  terms,  we  make  a  change  of  variables  so  that  the 
phase  of  Kl3  is  the  phase  of  a  delta  function.  Then,  we  determine  b,  by  the  criterion  that  in  order 
for  Kiti  to  best  approximate  a  delta  function,  its  amplitude  should  be  (27 r)~3. 


f  K(p,x)q(x)dx,  (3.29) 


The  term  A\,i. 


The  phase  of  Ki:i  is 

0i,i  (w,v,y’,p,x) 


u(Ti(y,y',x)  -  Ti (y,y’,p)) 

u[{\y  -x\  +  \x-  y'\)  -  (| y  -p\  +  \p-  y  |)] 
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(3.32) 


By  a  stationary-phase  calculation,  the  leading-order  contribution  to  K1 ;1  comes  from  the  stationary 
points  at  which  0  =  <9^01,1  =  <9y0i,i  =  <9y/0 i,i.  These  stationary  points  satisfy 

0  =  du</>ltl  =  [(Jy^*|  +  \^-y'\)  -  (\y-p\  +  \p-y\)] 

o  =  dy(j)  1,1  =  iy-^x)  -  (y-p) 

0  =  <9y/0M  =  (*  -  y')  -  (p  -  y')  (3.33) 

where  the  hats  denote  unit  vectors.  The  rst  equation  of  (3.33)  says  that  p  must  lie  on  the  same 
equal- travel-time  ellipse  as  x;  the  second  and  third  equations  say  that  the  directions  from  y  to  x 
and  p  must  be  the  same.  Clearly  the  only  point  satisifying  all  these  conditions  is  p  =  x. 

We  note  that  the  fonn  taken  by  the  critical  conditiions  depends  on  the  measurement  geometry. 
The  above  analysis  assumes  a  rather  unusual  situation,  namely  that  the  transmitters  y'  and  receivers 
y  are  spread  continuously  over  a  three-dimensional  region.  A  more  common  arrangement  is  for 
the  sources  and  receivers  to  form  a  planar  array  antenna;  in  this  case  the  integration  in  (3.31)  is 
only  over  the  two-dimensional  array,  and  the  differentiations  to  detennine  the  critical  points  are 
only  two-dimensional.  For  an  array  antenna,  the  second  line  of  (3.33)  is  replaced  by  (y  —  x)T  = 

{y  —  p)t  where  the  subscripts  T  denote  the  projection  onto  the  two-dimensional  array  plane  for 
the  receiver.  Similarly  the  third  line  of  (3.33)  is  replaced  by  a  projection  onto  the  array  plane  of  the 
transmitter.  We  note  that  the  two-dimensional  projection  of  a  unit  vector  determines  the  unit  vector 
in  the  case  we  have  here,  in  which  we  know  the  unit  vector  is  pointing  downwards.  (Physically 
the  fact  that  the  two-dimensional  projection  detennines  the  unit  vector  corresponds  to  the  fact  that 
an  array  antenna  can  produce  a  steered  beam.)  Thus  the  critical  equations  for  the  two-dimensional 
array  con  guration  imply  the  equations  for  the  three-dimensional  measurement  geometry.  For 
simplicity  of  notation,  we  write  simply  y  and  y '  for  the  positions  of  the  array  elements,  keeping  in 
mind  that  these  might  vary  over  only  a  two-dimensional  surface. 


The  term  Ad, 2. 


The  phase  of  K]2  is 

<t>i,2{u,y,y\p,x)  =  uj(T2(y,y',x)  -  T!(y,y',p)) 

=  u[(\y  -  z\  +  \z  -  x\  +  \x  -  y  |)  -  (\y  -  p\  +  | p-  y'\)] 

(3.34) 

The  stationary  points  satisfy 

0  =  <9^01,2  =  (\y-^z  \  +  |2j-a;|  +  \x  —  y'\)  -  (\y  —  p\  +  \p-y'\) 

0  =  dy(j)  1,2  =  (3/j^fO  - 

0  =  <9y/0i,2  =  {x  -  y')  —  {p  -  y')  (3.35) 


The  second  equation  of  (3.35)  says  that  p  must  lie  along  the  line  joining  2  and  y.  If  this  is  the 
case,  then  \y  —  z\  —  \y  —  p\  =  \z  —  p\.  The  rst  equation  of  (3.35)  then  becomes 

\z  —  x\  +  \x  —  y'\  =  \z  —  p\  +  \p  —  y'\  (3.36) 
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which  shows  that  p  must  lie  on  the  same  ellipsoid  as  x.  The  third  equation  of  (3.35)  speci  es  that 
the  p  must  be  in  the  same  direction  from  y'  as  x;  thus  p  =  x.  in  other  words,  all  conditions  of 
(3.35)  are  satis  ed  when  p  =  x  and  x  lies  directly  behind  2  when  viewed  from  the  position  y. 
This  situation  produces  an  image  p  of  the  point  x  in  the  correct  position,  but,  because  Ad  ,2  is  an 
off-diagonal  tenn,  the  strength  of  the  image  may  be  incorrect  at  such  a  location. 


The  term  Ad, 3. 


The  phase  of  K lj3  is 

(t>i,3{u,y,y  ,p,x)  =  u(T3(y,y',x)  -  Ti(y,y',p)) 

=  u[(\y  -x\  +  \x-  z\  +  \z-  y'\)  -  (\y-p\  +  | p-  y'\)] 

(3.37) 

Arguments  similar  to  those  for  Ad  ,2  show  that  the  only  critical  point  occurs  in  the  case  when  x  is 
directly  behind  2  as  seen  from  y';  in  this  case,  the  point  p  =  x  is  a  critical  point.  Again  the  image 
of  such  a  point  appears  in  the  correct  location,  but  its  strength  may  be  incorrect. 


The  term  Ad, 4. 


The  phase  of  Ad, 4  is 

0i,4(w,y,y' ,P,x)  =  u(Ti{y,y',x)  -  Ti(y,y',p)) 

=  u[(\y  -  z\  +  \z  -  x\  +  \x  -  z\  +  \z  -  y  |)  -  (| y  -p\  +  \p-  y  |)]  (3.38) 

The  stationary  points  satisfy 

0  =  00,01,4  =  {\y  ~  z\  +  \z  -  x\  +  I*  -  z\  +  \z  -  y'\)  -  (\y  -p\  +  \p  -  y'\) 

0  =  dy(j)  1,4  =  - 

0  =  dy'01,4  =  (z  -  y')  -(p-  y')  (3.39) 

The  last  two  conditions  of  (3.39)  say  that  p  must  he  on  the  lines  joining  y  to  2  and  y'  to  2.  Unless 
y  =  y',  this  implies  that  p  =  2;  the  rst  condition  of  (3.39)  can  then  be  satisi  ed  only  if  x  =  2. 
In  other  words,  this  tenn  can  produce  an  artifact  at  the  location  of  the  foreground  scatterer,  which 
is  not  a  location  in  which  we  are  interested. 


( y  -  p) 


The  term  K21. 


The  phase  of  Ad,i  is  the  phase  of  Ad, 2  with  y  and  y'  interchanged. 
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The  term  K22. 

The  phase  of  K2> 2  is 

4>2,2{u,v,v'  ,p,x)  =  u(r2{y,y',x)  -  r2{y,y',p )) 

=  u[(\y  -  z\  +  \z  -  x\  +  \x  -  y'\)  -  (| y  -  z\  +  \z  -  p\  +  \p  -  y'\)} 

(3.40) 

The  stationary  points  satisfy 

0  =  9^02,2  =  i\y  -  z\  +  \z  -  x\  +  \x  -  y'\)  -  (\y  -  z\  +  \z  -  p\  +  \p  -  y'\) 

0  =  dy(j) 2,2  =  (v^z)  - 

0  =  <9y/02,2  =  (*  -  y')  -  (p  -  y')  (3.41) 

The  second  condition  of  (3.41)  is  vacuous;  the  third  condition  says  that  p  lies  on  the  line  joining 
y ’  with  x.  The  rst  condition  (in  which  the  tenn  | y  —  z\  cancels)  says  that  p  must  lie  on  the  same 
ellipsoid  as  x.  The  only  such  point  is  p  =  x. 


iy-z) 


The  term  if2, 3. 


The  phase  of  if 2,3  is 

(t>2^(u,y,y  ,p,x)  =  uj(T3{y,y',x)  -  T2(y,y',p)) 

=  u[(\y  -  x\  +  \x  -  z\  +  \z  -  y  |)  -  (|  y  -  z\  +  \z  -  p\  +  \p  -  y\ )] 

(3.42) 


The  stationary  points  satisfy 

0  =  <9^02,3  =  {\yj^x\  +  \xj-z\  +  \z-y'\)-(\y-z\  +  \z-p\  +  \p-y'\) 

0  =  <9y02,3  =  (y^x)  -  (y^jz) 

0  =  dy>(j) 2,3  =  (z  -  y')  -  (p  -  y')  (3.43) 

The  second  equation  of  (3.43)  says  that  x  lies  directly  behind  2  when  viewed  from  y\  this  implies 
that  \x  —  y  \  —  \y  —  z\  =  \z  —  x\\  similarly  the  third  equation  of  (3.43)  says  thatp  lies  directly 
behind  z  when  viewed  from  y' ,  which  implies  that  \z  —  y' \  —  \p  —  y' \  =  —  \z  —  p\.  Thus  the  second 
and  third  equations  of  (3.43)  imply  that  the  rst  equation  reduces  to 

2\z  —  x\  —  2\z  —  p\.  (3.44) 

This  means  that  points  x  lying  along  the  line  joining  y  and  2  can  produce  artifacts  along  the  line 
joining  y'  and  2.  The  artifact  at  p  is  at  the  same  distance  from  2  as  the  scatterer  at  x. 
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The  term  K- 


The  phase  of  Ji2, 4  is 

02,4 (u,y,y',p,x)  =  u(r4(y,y\x)  -  T2(y,y',p)) 

=  u[(\y  -  z\  +  \z  -  x\  +  \x  -  z\  +  \z  -  y! |)  -  (| y  -  z\  +  \z  -  p\  +  \p  -  y'\)\ 

(3.45) 

The  stationary  points  satisfy 

0  =  cL02,4  =  i\y  -  z\  +  \z  -  x\  +  \x  -  z\  +  \z  -  y’\) 

___  -i]y_^z\  +  \z-p\  +  \p-y’\) 

0  =  <9y02,4  =  (VjZ^)  ~ 

0  =  <9y'02,4  =  (z  -  y')  -  (p  -  y')  (3.46) 

The  last  equation  of  (3.46)  says  that  p  lies  directly  behind  z  when  viewed  from  yl .  For  such  a 
point,  we  have  \z  —  y'\  —  \p  —  y'\  =  —  \z  —  p\  ,  which  converts  the  rst  equation  of  (3.46)  into 

2\x  —  z\  —  2\z  —  p\  (3.47) 

Thus  every  scatterer  x  produces  an  artifact  lying  directly  behind  2  when  seen  from  y'.  This  artifact 
is  at  the  same  distance  from  2  as  is  x. 


The  term  /y> , . 

The  phase  of  K-i4  is 

03,i {w,y,y',p,x)  =  u^y.y' ,x)  -  T3{y,y',p)) 

=  w[(\y-x\  +  | x-y'\)  -  (\y-p\)  -  {\p-z\  +  \z-y'\)\ 

(3.48) 

This  is  the  negative  of  the  phase  of  K \  3  with  x  and  p  interchanged.  Again  the  strength  of  scatterers 
directly  behind  2  can  be  incorrectly  reconstructed. 


The  term  /f3  2. 


The  phase  of  K:i  2  is 

(t>z,2{u,y,y'  ,p,x)  =  u(T2(y,y',x)-T3(y,y',p)) 

=  u[(\y  -  z\  +  \z  -  x\  +  \x  -  y'\) 

~(\y -p\)  ~  {\P~  z\  +  \z-y'\)\  (3.49) 

Again  this  is  the  negative  of  the  phase  of  iT2  3  with  x  and  p  interchanged;  this  tenn  can  cause 
scatters  lying  behind  2  from  y  to  appear  behind  2  when  viewed  from  y' . 
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The  term 


The  phase  of  K3  3  is 


03,3  (u,v,y',p,x)  =  u(T3{y,y',x)-T3(y,y',p)) 

=  u[(\y  -  x\  +  \x  -  z\  +  \z  -  y'\) 

-{\y-p\)-(\p-z\  +  \z-y'\)\  (3.50) 


The  stationary  points  satisfy 

0  =  cL03,3  =  (\y  -  x\  +  \x  -  z\  +  \z  -  y'\) 

-(I  y  -p[±_  lp  _  zl  + 12  _  y'\) 

0  =  <9y03,3  =  (y^x)  -  (y-j>) 

O  =  <9y'03,3  =  (z  -  y')  -  (z  -  y')  (3.51) 

In  the  rst  equation  of  (3.5 1),  the  tenn  \z  —  y'\  cancels;  we  see  that  p  must  lie  on  the  same  ellipsoid 
as  x  and  must  lie  in  the  same  direction  as  x  from  y.  The  only  such  point  is  p  =  x. 


The  term  A'3,4. 


The  phase  of  K:i  4  is 


03,4(^,?/!?/,!T'!a3)  =  w{n(y,y',x)-~T3(y,y',p)) 

=  u[(\y  -  z\  +  2\x  —  z\  +  \z  —  y'\) 

-( \y-p\)~(\p-z\  +  \z-y'\)\  (3.52) 

The  stationary  points  satisfy 

O  =  ct03,4  =  (\y-z\+2\x_-z\  +  \z-y’\)-(\y-p\  +  \p-z\  +  \z-y'\) 

0  =  <9y03,4  =  (VjZ^)  - 

0  =  <9y'03,4  =  0  -  y' )  -  (z  -  y')  (3.53) 

From  the  second  equation  of  (3.53),  we  see  that  y,  z,  and  p  all  lie  on  the  same  line,  which  implies 
that  \y  —  z\  —  \  y  —  p\  =  —\z  —  p\.  The  rst  equation  of  (3.53)  then  reduces  to2\x  —  z\  —2\p  —  z\. 

Thus  we  see  that  this  term  can  give  rise  to  artifacts  directly  behind  z  as  viewed  from  y. 


The  term  K44. 


04,4 {w,y,y',p,x)  =  u{Tt{y,y',x)  -  T4(y,y',p)) 

=  u[(\y  -  z\  +  2\z  -  x\  +  \z  -  y'\) 
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((I  y~z 


(3.54) 


The  stationary  points  satisfy 

0  =  <9^04,4 

0  b)y(f)  4,4 

0  —  dy'(f>  44 


+  2|z-p|  +  2  -t/  ))' 


((| y  -  z\  +  \z  -  x\  +  \x  -  z\  +  \z  -  y'\) 
-j\y^  zl  +j£_^p|  +  |p  -  z\  +  \z  -  y  |) 
(y^_z)  - 

(z  -  y')  -  (z  -  y>) 


(3.55) 


The  last  two  equations  of  (3.55)  are  vacuous;  the  rst  equation  reduces  to  \x  —  zj  =  \p  —  z\, 
which  says  merely  that  p  must  he  on  the  same  travel-time  sphere  about  2  as  does  x.  In  other 
words,  backprojection  of  a  path  in  which  both  the  incident  and  scattered  wave  bounce  off  of  2 
cannot  be  used  to  determine  the  location  of  a  scatterer  at  x.  This  is  because  the  scatterer  at  2  is 
isotropic:  after  a  wave  scatters  from  a  point  scatterer,  it  loses  all  information  about  the  direction 
from  which  it  came.  In  the  image,  this  term  causes  spherical  artifacts  around  2;  consequently  we 
omit  the  term  /i4  from  our  imaging  operator,  and  we  do  not  consider  tenns  of  the  fonn  if4j. 


Determination  of  the  bt 


We  have  found  that  the  imaging  operator  should  be  composed  of  three  terms: 

3  3  4  . 

m  =  E  D^fc](p)  ~  x)q{x)dx  (3.56) 

i= 1  i= 1  j= 1  J 

and  that  moreover,  if  we  restrict  our  attention  to  the  target  region  (avoiding  regions  behind  the 
scatterer  at  2),  only  the  diagonal  terms  Kht  contribute  (to  leading  order)  to  the  image.  We  have 
shown  above  that  this  imaging  operator  correctly  positions  scatterers  in  the  target  region. 

Next  we  turn  our  attention  to  the  scatterers’  strengths,  which  are  controlled  by  the  factors  6* 
appearing  in  Bt.  To  detennine  the  bi,  we  attempt  to  transfonn  each  Kt  l  into  a  delta  function. 
We  recall  that  a  delta  function  can  be  written  as  an  oscillatory  integral  in  the  form  5(p  —  x)  = 
(27 r)~3  J  exp | i (p  —  x)  ■  £]c/£.  Since  p  and  x  are  three-dimensional,  we  are  trying  to  express  Kr  l 
as  a  three-dimensional  integral.  This  means  that  our  measured  data  should  depend  at  least  three 
variables.  This  would  be  the  case,  for  example,  for  data  from  a  single  array  antenna,  where  the 
data  depends  on  t  and  two  array  coordinates  y.  If  more  data  is  available,  for  example  in  the  case 
in  which  we  have  a  separate  transmitting  and  receiving  arrays,  we  carry  out  the  analysis  below  for 
a  three-dimensional  subset  of  the  data  (say,  t  and  the  transmitter  coordinates  y)  and  then  simply 
integrate  over  the  remaining  variables. 

In  the  exponent  of  (3.3 1),  we  use  the  identity 

f(x)-f(p)=  I  ^/Gp  +  Kx  -p))d\  =  (x  -p)  ■  I  (V/)(p  +  \{x  -p))d\  (3.57) 
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to  write 


(3.58) 


u[Ti(y,y,x)  -  Ti{y,y',p)]  =  {x  -  p)  ■  Z\p,x,y,y',u ); 
explicitly,  the  E'  are  given  by 


liP,x,y,y',uj)  =  u  /  X7x'Ti(y,y\x') 


x,=p-\-X(x—p) 


dX 


When  p  =  x,  we  have 


=  ulp^y  +  p-^ 
£2{p,p,y,y',u)  =  uj[p^y  +  p^~z\ 
Z3{p,p,y,y',u)  =  u[p^z  +  P  -  y'} 


(3.59) 


(3.60) 


In  the  integral  (3.31)  for  and  K2, 2,  we  make  the  change  of  variables 

(u,y)  ->  £*  =  Er(p,p,y,y',Lv)]  (3.61) 

in  76 3  3  we  make  the  change  of  variables 

(w,y')  -»•  C  =  Z3(p,p,y,y',uy,  (3.62) 

This  transfonns  the  expressiion  (3.31)  for  A',  1  and  K2j 2  into 


Ki,%{p,x) 


d(uj,y)\ 
dC  ) 


(p,y,y',u) 


dCdy'  + 


(3.63) 


where  y  =  y(£)  and  cu  =  <u(£)  and  Et  denotes  a  smooth  error  tenn;  K.i:i  is  transformed  into  a 
similar  expression  except  that  the  integral  is  over  y  instead  of  y' . 


Equation  (3.63)  exhibits  the  point  spread  function  K  as  the  kernel  of  a  pseudodifferential  op¬ 
erator.  Pseudodifferential  operators  have  the  pseudolocal  property  [44],  i.e.,  they  do  not  move 
singularities  or  change  their  orientation.  It  is  immediately  clear  from  (3.63)  that  provided  the  Jaco¬ 
bian  | d(u,  y)/d£l\  is  nonzero,  the  leading  order  contribution  to  the  image  comes  from  the  points 
p  =  x. 


We  see  from  (3.63)  that  the  backprojection  weighting  function  b  should  be  chosen  as 

,  (a&)  Xi(p,y,y'x) 

bi(uj,p,y,y)  =  - 


(2vr)  ai(u,y,y',p) 


(3.64) 


where  \i  is  a  smooth  cutoff  function  that  prevents  division  by  zero  in  (3.64)  and  that  is  chosen  so 
that 


xi[p,y(£),y'M£)]  +  X2\p,v{€),vf  ,u{£)\  +  X3[p,y,y'(£)MO]  = 


(2*0 


(3.65) 
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in  as  large  a  region  of  £-space  (i.e.,  (a;,  r/)-space)  as  possible. 

The  Jacobian  detenninants  \d(u,  y)/d$,1 |  are  called  the  Beylkin  determinants  [?]  [8].  The  rst 
one  is 


d? 


=  det  p-y  +  p-  y' 


-ueP-~-  e1 
p-y 


-ueP-~-  e- 
p-y 


d(u,V) 

where  Pr  is  the  projection  operator  that  projects  a  vector  onto  the  plane  perpendicular  to  R: 

v  —  R(R  ■  v 


(3.66) 


Prv  = 


I  ill 


(3.67) 


and  where  e1  and  e2  denote  unit  vectors  tangent  to  the  receiving  array  antenna.  Similarly, 


_  ,  =  det  (  p  —  y  +  p  —  z  —l vP-rz-.e1  -ujP-rz-.e2 

d (lu,  y)  \  p~v  p~y 


of  ( _ _ 

777 - r  =  det  p  —  z  +  p  —  y'  — uP - -.  e 

d (u,  y)  \  p-y 


/i 


-uP — -.  e 
p-y' 


,/2 


where  e’1  and  e/z  denote  unit  vectors  tangent  to  the  transmitting  array. 


„/2 


(3.68) 

(3.69) 


The  determinants  of  (3.66),  (3.68),  and  (3.69)  can  be  calculated  easily.  These  determinants  are 
nonzero  because  their  column  vectors  are  linearly  independent:  for  example,  the  vector  p  —  y  + 
p  —  y'  points  from  the  antennas  towards  the  target,  whereas  Pjp^p  e'1  and  P^z^,  e/2  are  roughly 
tangent  to  the  antenna  array. 


Summary. 


For  the  case  of  a  single  point  scatterer  in  the  foreground,  the  imaging  operator  should  be 

3  3 

I(p)  =  Y.BA »‘c](p)  =  X  /  e'“lt-Tl{vvlv\(u,p,y,y’)g*(t,y,y’)dtdydy‘ 

3= 1  3= 1 J 

(3.70) 

where  the  r3  are  given  by  (3.27),  where  the  bj  are  chosen  as  in  (3.64),  and  where  the  Jacobian 
determinants  are  given  by  (3.66),  (3.68),  and  (3.69).  We  note  that  imaging  does  not  require  a  lot 
of  bookkeeping  in  the  sense  that  different  operators  do  not  need  to  be  applied  to  different  parts  of 
the  data.  Formation  of  the  imaging  operator  does,  however,  require  knowledge  of  the  foreground 
scatterer  and  does  require  that  the  backprojection  be  done  only  along  round-trip  paths  that  include 
a  direct  one-way  path  between  antenna  and  target. 
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3.4  Resolution 


With  the  bi  chosen  as  in  (3.64),  K  is  as  close  to  a  delta  function  as  possible  for  the  measurement 
geometry.  The  degree  to  which  it  approximates  a  delta  function  is  determined  by  the  support  of  x, 
which  is  in  turn  determined  through  (3.60)  by  the  overall  size  of  the  measurement  aperture. 

For  the  case  of  a  single  point  scatterer,  we  determine  as  follows  the  regions  in  ^-space  over 
which  we  have  data.  The  resolution  we  obtain  from  Kltl  is  detennined  by  the  region  in  Fourier 
space  covered  by  S1  as  y  and  y'  range  over  the  antennas  and  as  u  varies  over  the  bandwidth  of  the 
transmitted  waveform.  This  region,  which  is  determined  by  the  rst  equation  of  (3.60),  is  sketched 
in  Figure  3.2.  The  resolution  region  we  obtain  from  A'2,2  (sketched  in  Figure  3.3)  is  determined 
by  the  second  equation  of  (3.60),  and  the  resolution  region  obtained  from  AT3j3  (Figure  3.4)  is 
determined  by  the  third  equation  of  (3.60). 


y 


z 


Figure  3.2:  Region  in  Fourier  space  obtained  from  A\,i. 


3.5  Conclusions  from  the  theoretical  study 


We  have  exhibited  a  backprojection  imaging  method  that  makes  use  of  multipath  scattering  data 
from  point  scatterers  assumed  to  be  in  the  foreground  of  the  target.  We  nd  that  in  order  to  avoid 
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Figure  3.3:  Region  in  Fourier  space  obtained  from  Jl2, 2 

artifacts,  we  must  backproject  only  along  those  round-trip  paths  that  involve  a  direct  path  from  the 
target  to  the  antenna.  The  use  of  such  multipath  scattering  enhances  the  resolution  obtained  in  the 
image. 


3.6  Simulations 

Much  of  this  work  has  been  joint  work  with  Adam  Bojanczyk. 
The  simulation  work  proceeded  in  four  stages: 

1 .  Simulations  of  the  data  from  a  single  view 

2.  Reconstructions  from  simulated  single-view  data 

3.  Simulations  of  the  data  from  multiple  views 

4.  Reconstructions  from  simulated  multiple-view  data 
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Figure  3.4:  Region  in  Fourier  space  obtained  from  /v!  3 

3.6.1  Simulation  of  single-view  data 

For  the  simulations  we  used  a  frequency-domain  formula  for  multiple  scattering  from  point  scatter- 
ers  [1 1,  9].  The  code  allows  for  input  of  an  arbitrary  number  of  point  scatterers  located  at  arbitrary 
positions.  The  code  produces  stepped-frequency  data  for  a  frequency  range  that  can  be  input  by 
the  user.  This  code  produces  backscattered  data;  i.e.,  the  transmitter  and  receiver  are  assumed  to 
coincide. 


3.6.2  Reconstructions  from  a  single  view 

We  implemented  a  backprojection  algorithm,  which  is  e  xible  in  that  it  can  be  used  for  arbitrary, 
sparse  viewing  locations.  An  input  into  the  code  is  the  location  of  one  of  the  point  scatterers; 
knowledge  of  this  location  is  used  in  the  multipath  imaging  part  of  the  algorithm.  The  single¬ 
bounce  backprojection  from  this  known  scatterer  is  subtracted  out  of  the  image. 

The  algorithm  produces  a  backprojection  from  the  single-bounce  scattering  and  a  backprojec¬ 
tion  from  the  double-bounce  scattering.  Of  course,  since  the  algorithm  doesn’t  know  which  parts 
of  the  data  come  from  which  paths,  the  single -bounce  backprojection  contains  artifacts  due  to  the 
double-bounce  paths,  and  vice  versa. 

We  began  with  the  case  of  two  point  scatterers.  The  reconstructions  produced  in  this  case 
are  instructive.  Figure  2  shows  the  image  that  would  be  produced  by  a  standard  backprojection 
algorithm  that  does  not  account  for  multipath  scattering.  Recall  that  the  known  scatterer  has  been 
subtracted  out,  so  this  is  the  view  of  the  single  unknown  scatterer.  We  note  that  this  image  contains 
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input  data,  time  domain 


Figure  3.5:  A  time-domain  plot  of  (the  magnitude  of)  data  from  two  point  scatterers.  The  rst  two 
peaks  are  from  direct  scattering;  the  third  peak  is  from  the  wave  that  bounces  off  one  scatterer  onto 
the  second  and  then  propagates  back  to  the  radar 

no  cross-range  information  about  the  scatterer  location. 

Figure  3  shows  the  double-bounce  backprojected  image.  These  curves  are  ellipses  with  one 
focus  at  the  radar  and  the  other  focus  at  the  position  of  the  known  scatterer. 

Figure  4  shows  the  superposition  of  Figures  2  and  3.  In  this  gure,  we  can  nd  the  cross-range 
location  of  the  scatterer:  it  is  at  the  intersection  of  the  single-bounce  and  double-bounce  curves. 
We  note  also  artifacts  that  appear  behind  the  known  scatterer  (as  viewed  from  the  radar  position); 
this  is  predicted  by  the  theory. 

Backprojection  of  data  from  three  point  scatterers  is  much  more  complicated  and  is  dif  cult  to 
interpret.  For  this  reason  we  went  to  a  multiple-view  scenario. 


3.6.3  Simulation  of  multiple-view  data 

Simulation  of  multiple-view  data  was  done  by  including  rotation  of  the  point  scatterers.  A  range- 
angle  plot  of  simulated  data  is  shown  in  Figure  5.  We  see  sine  curves  from  each  scatterer,  together 
with  the  multiple-scattering  curve  on  top  (at  a  greater  range).  This  simulation  does  NOT  include 
attenuation  from  geometrical  spreading. 

We  have  run  numerical  experiments  with  the  correct  geometrical-spreading  attenuation  for 
three-dimensional  point  sources  and  scatterers;  the  amplitude  variation  then  makes  the  curves  more 
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dif  cult  to  see.  To  more  closely  simulate  the  experimental  data  from  the  range,  probably  the 
antennas  can  be  modeled  as  point  sources,  but  the  pipes  should  be  modeled  as  TWO-dimensional 
points,  so  that  the  geometrical  spreading  is  r~ 1/2  rather  than  r~l.  This  correction  has  not  been 
incorporated  into  the  code.  In  the  present  version  of  the  code,  attenuation  corrections  appropriate 
for  point  scatterers  in  three  dimensions  is  present  but  was  commented  out  to  produce  the  gures  in 
this  report. 


3.6.4  Multiple-view  reconstruction 

The  single-bounce  reconstruction  is  shown  in  Figure  3.10;  the  double-bounce  reconstruction  is 
shown  in  Figure  3.11,  and  the  sum  is  shown  in  F igure  3.12.  A  horizontal  slice  through  F igure  3.10 
is  compared  with  the  corresponding  slice  through  Figure  ??  in  Figure  3.13;  this  comparison  shows 
the  improvement  obtainable  by  including  the  double-bounce  backprojection. 

These  reconstructions  are  done  from  a  subset  of  the  data  shown  in  Figure  3.9,  which  does 
not  include  attenuation  from  geometrical  spreading.  Including  this  attenuation  is  problematic:  for 
well-separated  point  scatterers  in  three  dimensions,  the  returns  from  more  distant  scatterers  and 
the  multiple-scattering  returns  are  so  weak  that  multiple-scattering  artifacts  are  only  barely  visible 
in  the  single -bounce  reconstructions.  In  other  words,  multipath  artifacts  are  not  problematic  unless 
a)  the  scatterers  are  not  well-separated,  or  b)  the  three-dimensional  point  scatterer  assumption  is 
not  appropriate. 

When  attenuation  due  to  geometrical  spreading  is  included,  and  the  strength  of  the  double¬ 
bounce  reconstruction  is  boosted  suf  ciently  to  be  seen,  this  causes  problems.  In  particular,  boost¬ 
ing  the  strength  of  the  double -bounce  reconstruction  means  that  the  strength  of  the  artifacts  from 
single -bounce  scattering  are  also  automatically  boosted.  These  artifacts  tend  to  overshadow  the 
desired  image.  Persumably  these  dif  culties  would  be  alleviated  by  using  a  broader  bandwidth,  so 
that  the  data  more  closely  approximates  the  assumptions  of  the  theory. 


3.6.5  Work  with  data  from  the  radar  range 

The  rst  step  was  to  plot  the  data.  The  (magnitude  of)  the  data  is  shown  as  a  range-angle  plot  in 
Figures  3.14  (3  pipes),  3.15  (1  pipe),  and  3.16  (2  pipes).  Clearly  the  1-  and  2-pipe  data  is  much 
noisier  than  the  3-pipe  data.  All  three  datasets  show  scintillation,  the  interference  effects  from 
scatterers  in  the  same  or  nearby  range  cells.  This  effect  was  much  weaker  in  the  simulated  data 
above,  because  those  examples  are  from  scatterers  well-separated  in  range. 

The  reconstructed  image  from  the  3 -pipe  data,  using  single-bounce  backprojection  only,  is 
shown  in  Figure  8.  The  two  closer  scatterers  are  seen  clearly,  but  the  third  scatterer  is  obscured 
by  multipath  scattering  artifacts.  This  reconstruction  is  comparable  to  the  one  produced  by  Justin 
Bracken’s  code,  but  is  done  by  backprojection  rather  than  a  Fourier  transfonn. 
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A  double-bounce  backprojection  has  not  yet  been  successful,  because  I  have  not  yet  succeeded 
in  nding  the  location  of  one  of  the  scatterers  and  subtracting  it  out.  I  have  obtained  an  estimate  of 
the  location  of  of  the  strongest  scatterer,  but  subtracting  out  a  point  scatterer  at  that  location  does 
not  remove  it.  This  may  be  because  the  strongest  scatterer  was  a  pipe  suf  ciently  large  that  it  is 
not  behaving  as  a  point  scatterer. 


3.6.6  Issues  for  the  future 


The  following  issues  remain. 


1 .  Put  in  the  correct  geometrical  spreading  factors  into  the  data  simulator  and  into  the  recon¬ 
struction  codes. 

2.  Determine  whether  increasing  the  bandwidth  mitigates  the  single-bounce  artifacts  in  the 
double-bounce  reconstruction.  The  theory  predicts  that  for  suf  ciently  wide  bandwidth, 
singularities  (points  and  edges)  will  be  reconstructed  at  the  correct  positions,  and  the  rest  of 
the  image  (the  artifacts)  will  be  smoother.  How  wide  must  the  bandwidth  be  in  order  to  see 
these  effects? 

3.  Simulated  data  for  3  point  scatterers  does  not  look  like  the  experimental  range  data  for  3 
point  scatterers  -  why?  Is  this  only  because  good  estimates  for  the  pipe  positions  have  not 
been  put  into  the  code?  Or  are  the  waves  bouncing  more  than  once?  Is  there  a  resonance 
phenomenon  because  the  pipes  are  so  closely  spaced? 

4.  Write  code  for  nding  the  positions  of  the  closest  point  scatterers;  incorporate  this  into  the 
reconstruction  code. 

5.  Figure  out  how  to  subtract  the  strongest  scatterer  if  it  is  not  point-like.  It  may  be  possible  to 
treat  it  as  a  circular  superposition  of  point-like  scatterers,  nd  its  radius  using  an  optimization 
approach,  and  subtract  out  an  exact  solution  for  a  cylinder.  A  question  then  arises  of  its  effect 
on  the  multiple  scattering  processes;  will  the  multiple  scattering  from  such  a  large  object  be 
suf  ciently  similar  to  that  of  a  point  that  we  can  expect  the  theory  to  be  applicable? 

6.  Address  issues  of  noise  in  the  2-pipe  data. 

7.  Explore  alternative  approaches  for  dealing  with  multiple  scattering,  such  as  for  example 
using  the  exact  solution  for  multiple  scattering  from  point  sources  together  with  an  opti¬ 
mization  method  to  best  t  the  data. 
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Figure  3.6:  The  single-bounce  backprojected  image  of  a  point  scatterer,  with  the  true  locations 
shown  as  yellow  circles.  The  radar  position  is  the  red  cross  superimposed  on  the  scale  at  the 
bottom.  The  multipath  scattering  appears  as  a  second  scatterer  located  farther  away  from  the  radar. 
Note  that  no  cross-range  information  can  be  obtained  from  this  single  view. 
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Figure  3.7:  The  double-bounce  backprojected  image  from  two  point  scatterers,  with  the  true  loca¬ 
tions  shown  as  yellow  circles.  In  this  image,  the  single-bounce  scattering  causes  an  artifact. 
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Figure  3.8:  The  backprojected  image  from  two  point  scatterers;  this  is  a  superposition  of  Figures 
2  and  3.  The  true  locations  are  shown  as  yellow  circles.  See  the  text  for  discussion. 
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Figure  3.9:  Range-angle  plot  of  multiple-view  data  for  two  scatterers.  Here  views  are  taken  every 
.5°  over  angles  from  0  to  90°. 


Figure  3.10:  Single-bounce  reconstruction  from  10  views  spanning  90°  (each  9°  apart).  The  radar 
starts  from  the  top  and  moves  to  the  left  side. 
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Figure  3.11:  Double-bounce  reconstruction  from  the  same  dataset  as  Figure  3.10. 
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Figure  3.12:  The  sum  of  Figures  3.10  and  3.1 1  . 
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slice  through  image 


Figure  3.13:  The  red  line  is  a  horizontal  slice  through  Figure  3.10  (the  single-bounce  backpro- 
jection);  the  green  line  is  the  same  slice  through  Figure  3.12  (the  sum  of  the  single-bounce  and 
double-bounce  backprojections).  This  shows  the  improvement  from  adding  the  double-bounce 
backprojection. 
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Figure  3.14:  Range-angle  plot  of  the  (magnitude  of)  the  experimental  data  from  3  pipes. 
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Figure  3.15:  Range-angle  plot  of  the  (magnitude  of)  the  experimental  data  from  a  single  pipe. 
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Figure  3.16:  Range-angle  plot  of  the  (magnitude  of)  the  experimental  data  from  two  pipes. 
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Figure  3.17:  The  single-bounce  backprojected  image  from  the  3-pipe  data.  The  radar  moves  over 
12°  at  the  top. 
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